﻿using System;
using System.Collections.Generic;
using System.Text;
using SoftwareConsulting.BI11.Studio.InternalServices;
using SoftwareConsulting.BI11.NumericalMethods.Interfaces;
using SoftwareConsulting.BI11.NumericalMethods.Iterations;

namespace SoftwareConsulting.BI11.NumericalMethods.Regression
{    
    /// <summary>
    /// Класс полинома
    /// </summary>
    /// <remarks>c[0] + c[1] * x + c[2] * x^2 + ...</remarks>
    public class PolynomialFunction : IOneVariableFunction
    {
        private double[] _coefficients;//коэффициенты полинома

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="coefficients">Коэффициенты полинома</param>
        public PolynomialFunction(params double[] coefficients)
        {
            _coefficients = coefficients;
        }
        
        /// <summary>
        /// Возвращает полином, полученный прибавлением констатного слагаемого к существующему полиному
        /// </summary>
        /// <param name="poly"></param>
        /// <param name="newTerm"></param>
        /// <returns></returns>
        public static PolynomialFunction operator +(PolynomialFunction poly, double newTerm)
        {
            double[] coef = new double[poly._coefficients.Length];
            poly._coefficients.CopyTo(coef, 0);
            coef[0] += newTerm;
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает полином, полученный сложением двух существующих полиномов
        /// </summary>        
        public static PolynomialFunction operator +(PolynomialFunction poly1, PolynomialFunction poly2)
        {
            int n = Math.Max(poly1.Degree, poly2.Degree) + 1;
            double[] coef = new double[n];
            for (int i = 0; i < n; i++)
                coef[i] = poly1[i] + poly2[i];
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает коэффициент полинома в указанной позиции
        /// </summary>
        /// <param name="coefficientIndex"></param>
        /// <returns></returns>
        public double this[int coefficientIndex]
        {
            get 
            { 
                return coefficientIndex < _coefficients.Length ? _coefficients[coefficientIndex] : 0; 
            }
        }
        
        /// <summary>
        /// Возвращает полином с меньшей на единицу степенью, чем текущий
        /// </summary>
        /// <param name="polynomialRoot"></param>
        /// <returns>Полином, поделенный на (X - polynomialRoot)</returns> 
        /// <remarks>Проверок на правильность указания корня не ведется</remarks>
        public PolynomialFunction Deflate(double polynomialRoot)
        {
            int n = this.Degree;
            double remainder = _coefficients[n];
            double[] coef = new double[n];
            for (int k = n - 1; k >= 0; k--)
            {
                coef[k] = remainder;
                remainder = remainder * polynomialRoot + _coefficients[k];
            }
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает степень полинома
        /// </summary>
        public int Degree
        {
            get 
            { 
                return _coefficients.Length - 1; 
            }
        }

        /// <summary>
        /// Возвращает производную полинома
        /// </summary>
        /// <returns></returns>
        public PolynomialFunction Derivative()
        {
            int n = this.Degree;
            if (n == 0)
                return new PolynomialFunction(new double[] { 0 });            
            double[] coef = new double[n];
            for (int i = 1; i <= n; i++)
                coef[i - 1] = _coefficients[i] * i;
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает полином, полученный делением на указанное число
        /// </summary>
        /// <param name="dividend">полином</param>
        /// <param name="divisor">делитель</param>
        /// <returns></returns>
        public static PolynomialFunction operator /(PolynomialFunction dividend, double divisor)
        {
            return dividend * (1 / divisor);
        }

        /// <summary>
        /// Возвращает полином, полученный делением на указанный полином
        /// </summary>
        /// <param name="poly">делимое</param>
        /// <param name="divisor">делитель</param>
        /// <returns></returns>
        public static PolynomialFunction operator /(PolynomialFunction dividend, PolynomialFunction divisor)
        {
            return dividend.GetResultOfDivision(divisor)[0];
        }

        /// <summary>
        /// Возвращает полином-частное и полином-остаток в результате деления на указанный полином
        /// </summary>
        /// <param name="divisor">делитель</param>
        /// <returns></returns>
        public PolynomialFunction[] GetResultOfDivision(PolynomialFunction divisor)
        {
            PolynomialFunction[] result = new PolynomialFunction[2];
            int m = this.Degree;
            int n = divisor.Degree;
            if (m < n)
            {
                result[0] = new PolynomialFunction(new double[] { 0 });
                result[1] = divisor;
                return result;
            }
            double[] quotient = new double[m - n + 1];
            double[] coef = new double[m + 1];
            for (int k = 0; k <= m; k++)
                coef[k] = _coefficients[k];
            double norm = 1 / divisor[n];
            for (int k = m - n; k >= 0; k--)
            {
                quotient[k] = coef[n + k] * norm;
                for (int j = n + k - 1; j >= k; j--)
                    coef[j] -= quotient[k] * divisor[j - k];
            }
            double[] remainder = new double[n];
            for (int k = 0; k < n; k++)
                remainder[k] = coef[k];
            result[0] = new PolynomialFunction(quotient);
            result[1] = new PolynomialFunction(remainder);
            return result;
        }

        /// <summary>
        /// Возвращает интеграл при 0 для X=0
        /// </summary>
        /// <returns></returns>
        public PolynomialFunction Integral()
        {
            return Integral(0);
        }

        /// <summary>
        /// Возвращает интеграл при указанном значении для X=0
        /// </summary>
        /// <param name="valueFor0">значение при X=0</param>
        /// <returns></returns>
        public PolynomialFunction Integral(double valueFor0)
        {
            int n = _coefficients.Length + 1;
            double[] coef = new double[n];
            coef[0] = valueFor0;
            for (int i = 1; i < n; i++)
                coef[i] = _coefficients[i - 1] / i;
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает полином, полученный умножением на указанное число
        /// </summary>
        /// <param name="poly">полином</param>
        /// <param name="scale">число</param>
        /// <returns></returns>
        public static PolynomialFunction operator *(PolynomialFunction poly, double scale)
        {
            int n = poly._coefficients.Length;
            double[] coef = new double[n];
            for (int i = 0; i < n; i++)
                coef[i] = poly._coefficients[i] * scale;
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает полином, полученный перемножением двух полиномов
        /// </summary>
        /// <param name="poly1"></param>
        /// <param name="poly2"></param>
        /// <returns></returns>
        public static PolynomialFunction operator *(PolynomialFunction poly1, PolynomialFunction poly2)
        {
            int n = poly1.Degree + poly2.Degree;
            double[] coef = new double[n + 1];
            for (int i = 0; i <= n; i++)
            {
                for (int k = 0; k <= i; k++)
                    coef[i] += poly2[k] * poly1[i - k];
            }
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает корни полинома
        /// </summary>
        /// <returns></returns>
        public double[] CalculateRoots()
        {
            return CalculateRoots(MathExtended.DefaultNumericalPrecision);
        }

        /// <summary>
        /// Возвращает корни полинома с указанной точностью
        /// </summary>
        /// <param name="desiredPrecision"></param>
        /// <returns></returns>
        public double[] CalculateRoots(double desiredPrecision)
        {
            PolynomialFunction dp = this.Derivative();
            double start = 0;
            Random random = new Random();
            while (Math.Abs(dp.GetValue(start)) < desiredPrecision)
                start = random.NextDouble();
            PolynomialFunction p = this;
            NewtonZeroFinder rootFinder = new NewtonZeroFinder(this, dp, start);
            rootFinder.SetDesiredPrecision(desiredPrecision);
            List<double> rootCollection = new List<double>(Degree);
            while (true)
            {
                rootFinder.Evaluate();
                if (!rootFinder.IsDesiredPrecisionAchieved)
                    break;
                double r = rootFinder.BestZeroApproximation;
                rootCollection.Add(r);
                p = p.Deflate(r);
                if (p.Degree == 0)
                    break;
                rootFinder.Function = p;
                try 
                {
                    rootFinder.SetDerivative(p.Derivative());
                }
                catch (ArgumentException ex) 
                {
                    Trace.StartBlock("ArgumentException killed");
                    Trace.Message("ArgumentException catched and silently killed in CalculateRoots! Exception text:");
                    Trace.Message(ex.ToString());
                    Trace.EndBlock();
                };
            }
            return rootCollection.ToArray();
        }

        /// <summary>
        /// Возвращает полином, полученный вычитанием указанного числа
        /// </summary>
        /// <param name="poly">полином</param>
        /// <param name="subtrahend">вычитаемое</param>
        /// <returns></returns>
        public static PolynomialFunction operator -(PolynomialFunction poly, double subtrahend)
        {
            return poly - subtrahend;
        }

        /// <summary>
        /// Возвращает разность двух полиномов
        /// </summary>
        /// <param name="minuend">уменьшаемое</param>
        /// <param name="substrahend">вычитаемое</param>
        /// <returns></returns>
        public static PolynomialFunction operator -(PolynomialFunction minuend, PolynomialFunction substrahend)
        {
            int n = Math.Max(minuend.Degree, substrahend.Degree) + 1;
            double[] coef = new double[n];
            for (int i = 0; i < n; i++)
                coef[i] = minuend[i] - substrahend[i];
            return new PolynomialFunction(coef);
        }

        /// <summary>
        /// Возвращает как строку
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            StringBuilder sb = new StringBuilder();
            bool firstNonZeroCoefficientPrinted = false;
            for (int n = 0; n < _coefficients.Length; n++)
            {
                if (_coefficients[n] != 0)
                {
                    if (firstNonZeroCoefficientPrinted)
                        sb.Append(_coefficients[n] > 0 ? " + " : " ");
                    else
                        firstNonZeroCoefficientPrinted = true;
                    if (n == 0 || _coefficients[n] != 1)
                        sb.Append(_coefficients[n]);
                    if (n > 0)
                        sb.Append(" X^" + n);
                }
            }
            return sb.ToString();
        }

        /// <summary>
        /// Возвращает значение полинома при указанном значении X
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public double GetValue(double x)
        {
            int n = _coefficients.Length;
            double result = _coefficients[--n];
            while (n > 0)
                result = result * x + _coefficients[--n];
            return result;
        }

        /// <summary>
        /// Возвращает значение полинома и его производной при указанном значении X
        /// </summary>
        /// <param name="x"></param>
        /// <returns>первый элемент массива - значение, второй элемент массива - производная</returns>
        public double[] CalculateValueAndDerivative(double x)
        {
            int n = _coefficients.Length;
            double[] result = new double[2];
            result[0] = _coefficients[--n];
            result[1] = 0;
            while (n > 0)
            {
                result[1] = result[1] * x + result[0];
                result[0] = result[0] * x + _coefficients[--n];
            }
            return result;
        }
    }
}
